#################################################################
## Purpose: REPLICATION for "Criminal Victimization and Agency Attitudes in Mexico"
## Author: Cassy Dorff, 2024
## Step 3: Calc pred prob + create all model figures (MS & Appendix)
#################################################################

#################################################################
rm(list=ls())
source(paste0("~/Desktop/replication/0_setup.R"))

# load the data
load(file=paste0(pathData,"1_nationalsample_data.rda"))
load(file=paste0(pathData,"2_oversample_data.rda"))

# load models data
load(file=paste0(pathData,"model1_results.rda")) 
load(file=paste0(pathData,"model2_results.rda")) 
load(file=paste0(pathData,"model3_results.rda")) 
load(file=paste0(pathData,"model4_results.rda")) 
load(file=paste0(pathData,"model5_results.rda")) 
load(file=paste0(pathData,"model6_results.rda")) 
load(file=paste0(pathData,"model7_results.rda")) 
load(file=paste0(pathData,"model8_results.rda")) 
load(file=paste0(pathData,"model9_results.rda")) 
load(file=paste0(pathData,"model10_results.rda"))

# short name
ndata <- national

# set ref to largest group (vote)
ndata$ref<-relevel(ndata$dv2, ref="Vote")
over$ref<-relevel(over$dv2, ref="Vote")

# data for net
no_direct_vic <- ndata %>%
  filter(victim_score_combined == "0")

no_direct_vic_over <- over %>%
  filter(victim_score_combined == "0")
#################################################################

#################################################################
# calculate predicted probabilities from 
# model 1 using a counterfactual data.frame
# vary victimization at each level (0,1,2)
# average by victimization variable

pred_model1 <- predictions(
  model1,
  type = "probs",
  by = "victim_score_combined",
  newdata = datagridcf(victim_score_combined = 0:2))

# plot predicted probabilities 
model1_pred_plot <- 
  ggplot(pred_model1, 
         aes(x = victim_score_combined, y = estimate)) +
  geom_point() +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  xlab("Direct Victimization") + ylab("Predicted Probability") +
  scale_x_discrete(labels=c("0" = "No Direct \n Victimization", "1" = "Low-level \nVictimization",
                                "2" = "High-level \nVictimization")) +
  facet_wrap(~factor(group, levels = c(
    "Don't Know", "Nothing is possible", 
    "Join self-defense groups",
    "Protest", "Vote"))) +
  theme_minimal() 

ggsave(filename = paste0(pathGraphics, "model1_pred_plot.pdf"), plot = model1_pred_plot, dpi=350, width = 8, height=5)
#################################################################

#################################################################
# calculate predicted probabilities from 
# model 2 using a counterfactual data.frame
# vary victimization at each level (0,1,2)
# average by victimization variable

pred_model2 <- predictions(
  model2,
  type = "probs",
  by = "victim_score_combined",
  newdata = datagridcf(victim_score_combined = 0:2))

# plot preds
pred_model2_plot <- 
  ggplot(pred_model2, 
         aes(x = victim_score_combined, y = estimate)) +
  geom_point() +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  xlab("Direct Victimization") + ylab("Predicted Probability") +
  scale_x_discrete(labels=c("0" = "No Direct \n Victimization", "1" = "Low-level \nVictimization",
                            "2" = "High-level \nVictimization")) +
  facet_wrap(~factor(group, levels = c(
    "Don't Know", "Nothing is possible", 
    "Join self-defense groups",
    "Protest", "Vote"))) +
  theme_minimal() 

ggsave(filename = paste0(pathGraphics, "model2_pred_plot.pdf"), plot = pred_model2_plot, dpi=350, width = 8, height=5)
#################################################################

#################################################################
# calculate predicted probabilities from 
# model 3 using a counterfactual data.frame
# vary network victimization at each level (0,1)
# average by net victimization variable

pred_model3 <- predictions(
  model3,
  type = "probs",
  by = "victim_net_only_bin",
  newdata = datagridcf(victim_net_only_bin = 0:1))

# plot predicted probabilities 
model3_pred_plot <- 
  ggplot(pred_model3, 
         aes(x = victim_net_only_bin, y = estimate)) +
  geom_point() +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  xlab("Network Victimization") + ylab("Predicted Probability") +
  scale_x_discrete(labels=c("0" = "No \n Victimization", "1" = "Network \nVictimization")) +
  facet_wrap(~factor(group, levels = c(
    "Don't Know", "Nothing is possible", 
    "Join self-defense groups",
    "Protest", "Vote"))) +
  theme_minimal() 

ggsave(filename = paste0(pathGraphics, "model3_pred_plot.pdf"), plot = model3_pred_plot, dpi=350, width = 8, height=5)
#################################################################

#################################################################
# calculate predicted probabilities from 
# model 4 using a counterfactual data.frame (with more controls)
# vary network victimization at each level (0,1)
# average by net victimization variable

pred_model4 <- predictions(
  model4,
  type = "probs",
  by = "victim_net_only_bin",
  newdata = datagridcf(victim_net_only_bin = 0:1))

# plot predicted probabilities 
model4_pred_plot <- 
  ggplot(pred_model4, 
         aes(x = victim_net_only_bin, y = estimate)) +
  geom_point() +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  xlab("Network Victimization") + ylab("Predicted Probability") +
  scale_x_discrete(labels=c("0" = "No \n Victimization", "1" = "Network \nVictimization")) +
  facet_wrap(~factor(group, levels = c(
    "Don't Know", "Nothing is possible", 
    "Join self-defense groups",
    "Protest", "Vote"))) +
  theme_minimal() 

ggsave(filename = paste0(pathGraphics, "model4_pred_plot.pdf"), plot = model4_pred_plot, dpi=350, width = 8, height=5)
#################################################################

#################################################################
# calculate predicted probabilities from 
# model 5 using a counterfactual data.frame
# vary victimization at each level (0,1,2)
# average by victimization variable
# oversample

pred_model5 <- predictions(
  model5,
  type = "probs",
  by = "victim_score_combined",
  newdata = datagridcf(victim_score_combined = 0:2))

# plot preds
pred_model5_plot <- ggplot(pred_model5, 
                           aes(x = victim_score_combined, y = estimate)) +
  geom_point() +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  xlab("Direct Victimization") + ylab("Predicted Probability") +
  scale_x_discrete(labels=c("0" = "No \n Victimization", "1" = "Low-level \nVictimization",
                            "2" = "High-level \nVictimization")) +
  facet_wrap(~factor(group, levels = c(
    "Don't Know", "Nothing is possible", 
    "Join self-defense groups",
    "Protest", "Vote"))) +
  theme_minimal()

ggsave(filename = paste0(pathGraphics, "model5_pred_plot.pdf"), plot = pred_model5_plot, dpi=350, width = 8, height=5)
#################################################################

#################################################################
# calculate predicted probabilities from 
# model 6 using a counterfactual data.frame
# vary network victimization at each level (0,1)
# average by network victimization variable
# oversample

pred_model6 <- predictions(
  model6,
  type = "probs",
  by = "victim_net_only_bin",
  newdata = datagridcf(victim_net_only_bin = 0:1))

# plot preds
pred_model6_plot <- 
  ggplot(pred_model6, 
         aes(x = victim_net_only_bin, y = estimate)) +
  geom_point() +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  xlab("Network Victimization") + ylab("Predicted Probability") +
  scale_x_discrete(labels=c("0" = "No \n Victimization", "1" = "Network \nVictimization")) +
  facet_wrap(~factor(group, levels = c(
    "Don't Know", "Nothing is possible", 
    "Join self-defense groups",
    "Protest", "Vote"))) +
  theme_minimal()

ggsave(filename = paste0(pathGraphics, "model6_pred_plot.pdf"), plot = pred_model6_plot, dpi=350, width = 8, height=5)
#################################################################

#################################################################
# ROBUSTNESS
# model 7 // more controls

pred_model7 <- predictions(
  model7,
  type = "probs",
  by = "victim_score_combined",
  newdata = datagridcf(victim_score_combined = 0:2))

# plot preds
pred_model7_plot <- ggplot(pred_model7, 
                           aes(x = victim_score_combined, y = estimate)) +
  geom_point() +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  xlab("Direct Victimization") + ylab("Predicted Probability") +
  scale_x_discrete(labels=c("0" = "No Direct \n Victimization", "1" = "Low-level \nVictimization",
                            "2" = "High-level \nVictimization")) +
  facet_wrap(~factor(group, levels = c(
    "Don't Know", "Nothing is possible", 
    "Join self-defense groups",
    "Protest", "Vote"))) +
  theme_minimal()

ggsave(filename = paste0(pathGraphics, "model7_pred_plot.pdf"), plot = pred_model7_plot, dpi=350, width = 8, height=5)
#################################################################

#################################################################
# combine main models into one data.frame to plot comparison

# all preds all models
pred_model1$Model <- "model1" #h1
pred_model2$Model <- "model2" #h1 + controls (ms)

model1_model2_pred <- bind_rows(pred_model1, pred_model2)

model1_model2_plot <- ggplot(model1_model2_pred, 
                          aes(x = victim_score_combined, y = estimate)) +
                          geom_pointrange(
                            aes(ymin = conf.low, ymax = conf.high, colour=Model),
                            position = position_dodge(width = .3 )) +
                          xlab("Direct Victimization") + ylab("Predicted Probability") +
                          scale_x_discrete(labels=c("0" = "No Direct \n Victimization", 
                                                    "1" = "Low-level \nVictimization",
                                                    "2" = "High-level \nVictimization")) +
                          scale_color_manual(labels = c("Model 1", "Model 2"),
                                                    values = c("darkgrey", "black")) +
                          facet_wrap(~factor(group, levels = c(
                                      "Don't Know", "Nothing is possible", 
                                      "Join self-defense groups",
                                      "Protest", "Vote"))) +
                          theme_minimal() +
                          theme(legend.position="bottom")

ggsave(filename = paste0(pathGraphics, "model1_model2_pred_plot.pdf"), plot = model1_model2_plot, dpi=350, width = 8, height=5)
#################################################################

#########################################
# combine network models

# all preds both network models
pred_model3$Model <- "model3" #h2 net
pred_model4$Model <- "model4" #h2 + controls (ms)
model3_model4_pred <- bind_rows(pred_model3, pred_model4)

model3_model4_plot <- ggplot(model3_model4_pred, 
                              aes(x = victim_net_only_bin, y = estimate)) +
                              geom_pointrange(aes(ymin = conf.low, ymax = conf.high, colour=Model),
                                position = position_dodge(width = .3 )) +
                              xlab("Network Victimization") + ylab("Predicted Probability") +
                      scale_x_discrete(labels=c("0" = "No \n Victimization", "1" = "Network \nVictimization")) +
                      scale_color_manual(labels = c("Model 3", "Model 4"),
                      values = c("darkgrey", "black")) +
                      facet_wrap(~factor(group, levels = c(
                                "Don't Know", "Nothing is possible", 
                                "Join self-defense groups",
                                "Protest", "Vote"))) +
                      theme_minimal() +
                      theme(legend.position="bottom")

ggsave(filename = paste0(pathGraphics, "model3_model4_pred_plot.pdf"), plot = model3_model4_plot, dpi=350, width = 8, height=5)

#################################################################
# all controls, main model in MS in red for appendix
# get predictions
pred_model8 <- predictions(
  model8,
  type = "probs",
  by = "victim_score_combined",
  newdata = datagridcf(victim_score_combined = 0:2))

pred_model9 <- predictions(
  model9,
  type = "probs",
  by = "victim_score_combined",
  newdata = datagridcf(victim_score_combined = 0:2))

pred_model10 <- predictions(
  model10,
  type = "probs",
  by = "victim_score_combined",
  newdata = datagridcf(victim_score_combined = 0:2))

# labels for plotting
pred_model7$Model <- "model7" #h1 + controls (app1)
pred_model8$Model <- "model8" #h1 + controls (app2)
pred_model9$Model <- "model9" #h1 + controls (app3)
pred_model10$Model <-"model10" #h1 + controls (app4)

# data for plotting
all_h1_models_pred <- bind_rows(pred_model2, pred_model7, pred_model8,  pred_model9, pred_model10)

model_h1_app_plot <- ggplot(all_h1_models_pred, 
                            aes(x = victim_score_combined, y = estimate)) +
  geom_pointrange(
    aes(ymin = conf.low, ymax = conf.high, colour=Model),
    position = position_dodge(width = .4 )) +
  xlab("Direct Victimization") + ylab("Predicted Probability") +
  scale_x_discrete(labels=c("0" = "No Direct \n Victimization", 
                            "1" = "Low-level \nVictimization",
                            "2" = "High-level \nVictimization")) +
  scale_color_manual(labels = c("Model 2", "Model 7", "Model 8",
                                "Model 9", "Model 10"),
                     values = c("tomato", "lightgrey", "grey", "darkgrey","dimgrey", "black")) +
  facet_wrap(~factor(group, levels = c(
    "Don't Know", "Nothing is possible", 
    "Join self-defense groups",
    "Protest", "Vote"))) +
  theme_minimal() +
  theme(legend.position="bottom")

ggsave(filename = paste0(pathGraphics, "appendix_h1_all_plot.pdf"), plot = model_h1_app_plot, dpi=350, width = 8, height=5)



